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Abstract 

We apply second order finite difference to calculate the lowest eigenvalues of the 
Helmholtz equation, for complicated non-tensor domains in the plane, using different 
grids which sample exactly the border of the domain. We show that the results obtained 
applying Richardson and Pade-Richardson extrapolation to a set of finite difference 
eigenvalues corresponding to different grids allows to obtain extremely precise values. 
When possible we have assessed the precision of our extrapolations comparing them 
with the highly precise results obtained using the method of particular solutions. Our 
empirical findings suggest an asymptotic nature of the FD series. In all the cases 
studied, we are able to report numerical results which are more precise than those 
available in the literature. 


1 Introduction 

Among the different methods for estimating the eigenvalues and eigenfunctions of the Lapla- 
cian on a finite region of the plane, finite differences (FD) is the simplest, although the ac- 
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curacy of the results obtained with this method is limited. In particular, for domains with 
reentrant corners with an angle of 7r/a, it is well known that the error of the FD eigenvalues 
is dominated by a behavior for h —)• 0 (h is the grid spacing). 

The so-called L-shaped membrane [a = 4/3] is a famous example which was studied 
long time ago by Fox, Henrici and Moler M- Because of the quite slow convergence of 
FD in this case {AE ^ those authors applied an alternative method, the method of 

particular solutions (MPS), and, exploiting all the symmetries of the problem, they were able 
to obtain the hrst 8 digits of the lowest eigenvalue of the L-shape correctly, Ei Ri 9.6397238. 
Interestingly, the paper also mentions a precise (unpublished) value obtained by Moler 
and Forsythe, Ei 9.639724, extrapolating the FD values obtained with very fine grids. 
Unfortunately, the extrapolation is neither named nor explained. 

A valuable discussion of the Richardson extrapolation of FD results for the eigenvalues 
of the Laplacian on two dimensional regions of the plane is contained in |21j . where it is 
pointed out that the correct exponents of the asymptotic behavior of Ei for h ^ 0 must be 
used in the extrapolation. 

The purpose of the present paper is to show that is it possible to obtain quite precise 
approximations to the eigenvalues of the Laplacian on a certain class of two dimensional 
domains (specihcally domains whose borders are sampled by the grid) by Richardson ex¬ 
trapolation of the FD results, provided that the asymptotic behavior of the FD eigenvalues 
for h —>■ 0 is taken into account correctly. 

The paper is organized as follows: in section[2]we provide a general discussion of Richard¬ 
son extrapolation, and its relation to the “method of deferred corrections”; in section [H we 
describe the practical implementation of the Richardson extrapolation used in this paper; 
in section |T] we present the numerical results obtained for different domains, comparing 
them with the best results available in the literature; finally, in section [5] we summarize our 
findings and discuss possible directions of future work. 


2 Richardson Extrapolation 

Richardson Extrapolation is interpolation of samples of a sequence Sn by a continuous 
function of a continuous variable z followed by extrapolation to z = 0 to approximate the 
limit of the sequence. The slowly convergent series for example, can be summed 

by taking the sequence of partial sums. Si, = fo be samples of a function in 

z = l/iy. In our application, the sequence is that of approximations to an eigenvalue by 
finite difference calculations whose asymptotic error is a series in some power of the grid 
spacing h; here z = hS [usually] or z = [for one singular application.] 

The history including many independent discoveries is reviewed by Brezinski [7] , Marchuk 
and Shaidurov [H], Sidi [33], Walz [37] and Joyce [30]. Christian Huyghens applied Richard¬ 
son Extrapolation to estimate tt to 35 decimals from the perimeters of a sequence of polygons 
with more and more sides inscribed in the unit circle. Richardson’s (1927) paper |29j con¬ 
tained a plethora of examples that was the Hrst comprehensive display of the power of 
extrapolation; he claimed no novelty but credited others including an obscure Russian lan¬ 
guage paper by Bogolouboff and N. Krylov 0 Richardson Extrapolation of eigenvalues is 
discussed in Pryce’s book on numerical solution of Sturm-Liouville problems m- 

Richardson Extrapolation has four steps. First, compute samples {/(h„)} of the function 
being extrapolated. Second, choose a set of basis functions {4>j{x)} - usually polynomials - 

^N. Bogolouboff and N. Krylov, On the Rayleigh’s principle in the theory of he differential equations 
of the mathematical physics and upon the Euler’s method in the calculus of variations, Acad, des Sci. de 
rUkraine, Classe, Phys. Math., tonne 3, fasc. 3 (1926). 
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for an approximation 


N 

fN{h)='^ aj4>j{h) ( 1 ) 

i=i 

The coefficients aj can always be computed by solving a matrix problem at a cost of 0{N^ ) 
operations, and this is necessary when the (pj are a mixture of polynomials and polynomials 
multiplied by powers of log(x), for example. However, it is faster to use Neville-Aitken 
interpolation to compute a two-dimensional array (’’Richardson Table”) of approximations 
of different N formed from different subsets of the full sample set {/(/in)}- This is cheaper 
than matrix-solving [0{N'^) floating point operations] though this is only a small virtue 
because of the speed of modern laptops. More important, extrapolation is credible only 
if its answers are independent of numerical choices such as N and subsets of the full set 
of samples. More precisely, a numerical answer is believable if and only if several different 
values of the numerical parameters yield the same answer to within the user chosen tolerance. 
The Richardson Table allows a quick search for such stable approximations. We shall return 
to this in analyzing each numerical example. 

Various conventions are employed. A popular one is to arrange the table as a lower 
triangular matrix with N samples of f{z), the function being approximated, as the first 
column: 


Rj,i = fizj) 


( 2 ) 


The simple recursion is 


Rj^k — 


{z Zj-k-i)Rj,k-i {z Zj)Rj-i^k-i 
fc+i 


k=j,{j + l),...N, j = l,2,...7V (3) 


Each entry in column A: is a polynomial of degree (fc — 1) which interpolates a subset of k 
samples. The basic step combines two polynomials that interpolate (fc — 1) points each to 
generate a polynomial that interpolates at the k points {zj-k+iT ■ ■ Both generators 
interpolate at the (fc — 2) points {zj-k+i, ■ ■ ■ Zj}, but only Rj^k-i interpolates at Zj while 
Rj-i^k-i does not, but interpolates at Zj-k-i-i- It is easy to verify that 


Rj,k{z — Zj — k+l) 


{Zj — k+l Zj — k—l)Rj^k—l {Zj — k-\-l Zj)Rj — i^k — l 


fizj-k+i) [using = Zj-k+i) = f{zj-k+i)] 


(4) 

(5) 

( 6 ) 


Rjjki^Z — Zj^ — 


Zj-k-l)Rj,k-l (.Zj Zj)Rj-i^k- 


(^j — k-\-l 

^j,k — l 


i^j Zj-k-l 
^3 ^j — k+1 

f(zj [using Rj,k-i(z = Zj) = f(zj)] 


(7) 

( 8 ) 
(9) 
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—fc+l) — 


i^j-k+l ^j-k 
- ^ 3 ) 


; 1 ) Rj ^k — 1 ( ^j — fc + 1 ^j ) Rj — 1, fc — 1 

•y - _ -y . , , , 


^j — k-\-l 


■Rj-i^k-1 


^j,k{^ — ^m) — 


^3 ^j — k-\-l 

f{Zj-k+l) 

{^m ^j — k—l)Rj,k—l {^m ^j^^j — l,k — l 

- - - -^^-—— -^-, 171 = J 

^3 ^j — k-\-l 

{Zm — Zj-k-l) — {Zm — Zj) ^ 

/ \Zm) 

Z] -z, ■ - 

~ ^j — k—1 


( 10 ) 

( 11 ) 

( 12 ) 

= j -k + 2,...j -I 




^3 ^j — k-\-l 

= f{Zm) 


(13) 

(14) 

(15) 


where we used Rj,k-i{zm) = Rj-i,k-i[zm) = /(-Zm) in the last lines. 

For Richardson Extrapolation, we set z = 0 and the table of polynomials becomes a 
lower triangular matrix of numbers. 

When z = 1/n, a reciprocal integer, Salzer gave a nice closed-form extrapolation formula 
in 1954 [30] as well as tables of the weights assigned to each sample in the final answer. 

Sidi gives some convergence proofs in Chapter 3 of his book [33] . It is known that 
Richardson Extrapolation is often exponentially (geometrically) convergent with the error 
of the diagonals and bottom rows of the table falling as exp(—gn) for some positive constant 
q even when the power series being extrapolated is factorially divergent, as usually true when 
the samples are of the trapezoidal rule for different grid spacings h and the associated series 
in powers oi z = h3 is the Euler-Maclaurin formula. A comprehensive theory is still lacking, 
however. 

Richardson Extrapolation is closely related to the “method of deferred corrections”, al¬ 
ternatively labelled “correction by higher order differences” in the (1983) book by Marchuk 
and Shaidurov [^. “Deferred corrections” also solves matrix problems that are the low 
order, usually second-order, discretization of the problem. Deferred corrections also pro¬ 
motes this low order approximation into a very high order approximation. In contrast to 
Richardson Extrapolation, which solves the low order problem repeatedly on a variety of 
different grids, deferred corrections uses only a single grid, and applies an iteration precon¬ 
ditioned by the low order discretization nail]. The residual is evaluated by a high order 
method; the accuracy of the converged iterative solution is equally high. One grid, instead 
of many, is obviously a significant advantage for deferred correction. The method can be 
applied to eigenvalue problems [36119] .This approach has become the standard way of gen¬ 
erating very high order time marching schemes to pair with spectral spatial discretizations. 
Dutt, Greengard and Rokhlin write, “We begin by converting the original ODE into the 
corresponding Picard equation and apply a deferred correction procedure in the integral for¬ 
mulation, driven by either the explicit or the implicit Euler marching scheme. The approach 
results in algorithms of essentially arbitrary order accuracy for both non-stiff and stiff prob¬ 
lems” |12] . Further developments of Picard integral/deferred correction time-marching can 
be found in [nnnns]. 

High order evaluation on a line in one dimension (time) is easy, but evaluating the resid¬ 
ual of a partial differential equation by, say, twelfth order finite differences, is a bookkeeping 
nightmare. The programming and debugging escalate rapidly when the domain is geomet¬ 
rically complicated. Furthermore, corner singularities may make higher order evaluation of 
the residual impossible without heroic measures [B] . For all the success of deferred correction 
in other applications, for eigenproblems in domains with corners Richardson Extrapolation 
is clearly the better way. 
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3 Implementation of Richardson extrapolation 

Suppose that we have calculated a given eigenvalue of the Laplacian on a certain domain 
using finite differences for a number of grids, which all sample the border, and with de¬ 
creasing grid spacings, hi > h 2 >■•■> h^- Only when ft, ^ 0 is the exact eigenvalue of 
the associated problem in the continuum obtained, although the eigenvalues obtained for 
different (finite) grid spacing an asymptotic behavior, which depends on ft; for the grid 
we may typically expect 


E- 


(fe) 


= Co -b ^ . 


’ h ^ 


(16) 


where ai < a 2 < ■ ■ ■ < o-n ■ However, logarithms and more exotic functions have arisen 
in other problems. The exact values of these coefficients will depend on the particular 
properties of the domain studied: in fact, while integer values of a are associated with the 
discretization of the problem (a = 2,4,...), rational values of a may also appear when 
reentrant corners are present (as for the case of the L-shape where ai = 4/3). 

Using eq. (HU) for all grids, and with basis functions (j)j, one obtains a system of linear 
equations 


= co(/o + Ci(/i(fti) -b C24‘2{hi) H-b C7v-i(/Ar-i(fti) + . ■. 

/o') 

E{ = Co(/o + Ci4>i{h2)c24>2ih2) H-b CN-l4>N-l{h2) + ... 

e[^'’ = Co(j)o + Ci(/i(ftAr) -b C24>2(hN) ^ -b CN-l4>N-l(hN) + ■ • ■ 


(17) 


where the unknowns are the coefficients Cj (/ = 0,1,..., — 1). 

In matrix form these equations take the form 


where 



^ Co ^ 


( \ 

R 

Cl 

= 



\ CN -1 / 




</o 

(/l(ftl) 

(j) 2 {hi) . 

(j)N- 


</o 

4 >iih 2 ) 

4 > 2 {h 2 ) ■ 

4 >N- 

-1(^2) 

</o 

4 >i[hN) 

</ 2 (^Af) ■ 

4 >n- 

-i(^Ar) 


The solution to Eqs. (ESI) is obtained as 


/ Co \ 


( \ 

Cl 

= R-^ 


\ cat-i y 


V ) 


(18) 


(19) 


( 20 ) 


where the extrapolated value of cq will provide an estimate of the exact eigenvalue. 

Cramer’s rule can be used to obtain the coefficients Cj without inverting the matrix R; 
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in particular 


Co 


In our numerical examples 


El ... 

E 2 (t>i{h2) ... 0Ar-i(/i2) 

En (j)l{hN) ... (j)N-l{hN) 

(t>0 </'l(^l) (j)N-lihl) 

4>0 (t>l{h2) ... 4‘N-l{h2) 

4>0 4>N-l{hjs[) 


<pjiz) = 2 “^ 


( 21 ) 


( 22 ) 


where ao = 1 and the aj are a monotonically increasing sequence of positive constants. 

When we apply eq. dm) to the different grids, we are implicitly assuming that h > 
hi > > /iTv, where h is the radius of convergence of the series. However, in general 

h is unknown and it will only be estimated once the first few coefficients cj have been 
approximated. For this reason inaccurate results could be obtained if the spacing of one of 
the grids falls outside the radius of convergence of the asymptotic series. This is a common 
problem also of perturbative series, which are known to be divergent in many cases. 

To avoid this problem, we can extrapolate by Fade rational approximation 




Co 


E n 7 o 

7 = 1 Cj-h;, 


^ + Ef=idjh: 


(23) 


For integer exponents, aj and 13j, and N = M, the choice a at = /Sat, would correspond 
to a diagonal Fade. In a general case, with N ^ M and rational exponents, we assume 
aN = 13 m- 

Using the different grids (in this case we use iV + M + 1 grids) we obtain the system of 
linear equations 

i;(i) = co + cih'^^+--- + CNh'^^-dih^i^E^^^ - dMh{^E^^'^ 

i;(2) = Co + cih’^^ + ■ ■ ■ + CNh^^ - dihl^ - dnh^^ 


j^(N+M+l) ^ co + Cih‘^^\M+l + --- + CNh^N\M+l-dl 

which can be cast in matrix form as 


^ Co ^ 


Cl 


( \ 




Cn 

— 


di 




\ dM / 


where 

/ 1 




_h3MEii) \ 


R = 

1 


hr 



(25) 


1 ’ 1 ’ 

uai 

"'N+M+1 ■ 

^Af+M+l 
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(26) 




\ Am / 

or using Cramer’s rule once again. 


4 Numerical results 

To apply the extrapolation schemes described in the previous section we need to calculate 
accurately the FD eigenvalues for a series of grids. We consider different domains, with 
borders which can be sampled by a square grid and with different reentrant angles. 

4.1 L-shaped domain 

We consider the L-shaped region V, = {|x| < 1, \y\ < 1} — {0 < a; < 1, 0 < y < 1}, represented 
in Fig. [H Using finite differences and a five-points approximation to the Laplacian, the 
Helmholtz equation on U is solved with Dirichlet boundary conditions on 912 for a series of 
grids with an increasing number of points. We have exploited the symmetry of the domain, 
to obtain separately the even and odd modes of the L-shape. 

Our numerical calculations consist of two sets: 

• A calculation of the lowest eigenvalue of the L, using 124 grids with spacing h = I/Nq 
and Nq = 10,..., 133. The finite difference results of this set are obtained using the 
’’Conjugate Gradient Method” (CGM), as described in Ref. [5S], and they are accurate 
to 220 digits; 

• A calculation of the lowest 100 eigenvalues of the L, using 100 grids with spacing 
h = I/Nq and Nq = 10,..., 109. The finite difference results of this set are obtained 
using the internal Mathematica command Eigenvalues and they are accurate to 60 
digits. 

In Table [T] we report the available estimates of the lowest eigenvalue of the L-shape in 
the literature, including the results of the present work. 

As we have mentioned before, the convergence of the numerical results is affected by the 
presence of a reentrant corner and the finite-difference eigenvalue E{h) behaves for h ^ 0 

as[T0l[2T] 


E{h) = E{0) + -f ... 


(27) 


where E{0) is the eigenvalue of the Laplacian in the continuum. For the related problem of 
a H-shaped membrane, Donnelly m conjectured the asymptotic behavior 


E(h) = E{0) + -b bh^ + + dh'^ + ... 


(28) 


for the fundamental eigenvalue El- This behavior was also used by Christiansen and Petersen 
[8] to perform a Richardson extrapolation of the finite difference results for the L-shape (see 
Table [ID. 


^Since the H-shaped domain contains the same reentrant angle of the L-shape, we assume the same 
asymptotic law for both domains. 
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Table 1: Available estimates of the lowest eigenvalue of the L-shape (smaller fonts are used 
for the last three values, to allow fitting the results in the column). 



El 

Reid and Walsh [2S| 

9.63972 

Fox, Henrici and Moler [14] 

9.6397238 

Mason [25] 

9.6397 

Sideridis [32] 

9.6395 

Schiff [31] 

9.659 

Christiansen and Petersen [5] 

9.63972383991 

Still [35] 

9639723ff 

Betcke and Trefethen [^ 

9.6397238440219 

Amore [I] 

9.6397238440 

Yuan and He [1^ 

9.63972384^1 

this AiVOrk (R.ichardson) 

9.63972384402194105271145926236482315626728952582190645 

this work (Pade-R.ichardson) 

9.6397238440219410527114592623648231562672895258219064561095797005640 

this work (mps) 

9.639723844021941052711459262364823156267289525821906456109579700564036 


Figure 1: L-shaped region 


The results obtained extrapolating the FD sequences can be compared with the precise 
results obtained with the ’’method of multiple solutions” (MPS)[T4]. Table[2]reports the first 
25 eigenvalues of the L-shape obtained with the MPS (for the case of the first eigenvalue we 
have used 545 points evenly spaced, which allow one to obtain 70 digits of precision, for the 
remaining cases we have used 425 points, which allows an accuracy of about 50 digits). The 
eigenvalues marked with f are known exactly and correspond to modes of a square. The 
MPS has been implemented in Mathematica 10 [38], taking advantage of Mathematica’s 
ability to work with arbitrary precision numbers or with a large number of digits (in our 
case typically numbers are specified to 100 digits). 

We will use these values to establish the accuracy of the approximate values of 
obtained by applying four different extrapolation schemes, differing in the choice of the 
exponents: 











• Extrapolation i 


E{h) = £;(0) + ^cW/r^”RiE(0)+4*^/i2 + 4*^/r'‘ + 0(/r®) (29) 

n—1 

• Extrapolation ii 

OO 

E{h) = E(0) + ^c^")/i"a:E;(0) + cf^/i2 + 4")/i3 + 0(/i4) (30) 

n—1 

• Extrapolation iii (Donnelly, Ref. [10) 1 


OO 


E{h) 

= E(o)+^[4fhp"-"/=+4':’p" 

Yi — 


• Extrapolation iv 

« E{0) + + o{h^) 

(31) 

E{h) 

OO 

= E;(0) + ^cl“)/i2("+l)/3 

Th - 1 



« E;(0) + + cf + cf + 0(/ll°/3) 

(32) 


The first two schemes only use integer exponents and are expected to be accurate only 
for the modes of the L-shape which are also modes of the square. 

Figure [2] displays the error fQj- lowest eigenvalue of the L-shaped 

region, using the third and fourth extrapolation schemes. Here 

Aa = 124) ) _£,(MPS) I ^33^ 

Ab = _7^(fe-i,i24)(-^^)| (-34^ 

where the superscripts (iii) and (iv) refer to the series used and the ED eigenvalues are 
accurate to 220 digits. The values are the analogous of but using ED eigenvalues 
are accurate to 60 digits. 

The approximations obtained with the first two schemes, which do not use rational 
exponents, are very poor for this mode. 

In particular, the extrapolated values in the four cases are 

• Extrapolation i 

El w 9.639 8 (35) 

• Extrapolation ii 

El 9.6397 327 (36) 

• Extrapolation iiii 

El Ri 9.639723844021 1929465 (37) 

• Extrapolation iv (corresponding to the minimum in Fig. [2]) 

El ^ 9.639723844021941052711459262364823156267289525821906456 458 (38) 
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Remarkably, the fourth scheme provides the first 55 digits of Ei for the L-shape correctly, 
suggesting that the the exact asymptotic behavior of the finite difference eigenvalues, for 
h ^ 0, is E{h) = E{0) + ZZi 

In correspondence to the minimum of Fig. [2] we have calculated the first few coefficients 
of the asymptotic series for the eigenvalue of the fundamental mode; the expansion reads 
(underlined digits are expected to have converged) 

E(h) ss 9.63972384402194105271145926236482315626728952582190645 6 

+ 2.19759909080385142157537952672409583683648557094 5 

- 5.254349649878412271190008297029240841285038851 0 

- 0.045716100985365949827658978449794728350032 8 

- 1.946468144036811059220897747699440650587 7 

+ 1.125074754927755172836371946813777186 1 

- 0.2147544087374345021476728527871998 5 

+ 0.35588422353456505262712958896229 4 

+ 0.006403070910486707732478038349 7 h® 

+ 0.038286091425541761563564936 0 h^o/® 

- 0.0730523282127573068239088 6 h^^/® + ... (39) 

The behavior of the error in Fig. [2] suggests that the FD series is asymptotic. Therefore, 
if one picks a set of grids with spacings hi > /12 > ..., it is convenient to perform an 
extrapolation using the grids up to a given spacing hjv where the error reaches a minimum. 

This behavior, however, does not limit the number of accurate digits of the eigenvalue 
that one can obtain using the Richardson extrapolation. This is illustrated in Figs. [3] and H) 
the first figure is obtained extrapolating the FD results of a set with smallest spacing hmin 
and determining the minimum error over the extrapolated eigenvalue (which will correspond 
to the minimum observed in Fig. [21). In this case we observe that the number of accurate 
digits of the extrapolated eigenvalue grows linearly for A^o ^ 1- Of course this behavior will 
be lost when the number of digits of the FD eigenvalue is not sufficient (see for example, 
the last curve of Fig. [21 where the FD eigenvalue are only accurate to 60 digits). Fig. [4] 
illustrates the fact that, as hmin gets smaller and smaller, the number of grids used in the 
optimal extrapolation also grows linearly. 

In Fig. |S]we have applied the Fade-Richardson extrapolation to calculate the error over 
the fundamental eigenvalue of the L. Here F(^424) indicates the diagonal Fade with 2fc + 1 
coefficients, which uses the grids going from 124 —2fc to 124. The horizontal line corresponds 
to the lowest error obtained with the Richardson extrapolation, i.e. to the minimum of Fig.|21 
The errors are obtained using as a reference the precise estimate obtained using the MFS 
with 545 points distributed on the border, which is expected to have at least 70 correct 
digits (see Table [J). 

The result obtained with the Fade-Richardson extrapolation contains 13 extra digits of 
accuracy with respect to the result obtained with the Richardson extrapolation alone!! 

The same analysis can be carried out for the eigenvalue of the first excited mode of the 
L-shaped membrane, which is odd with respect to reflection about the line y = x] also in 
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Figure 2: Error over the first eigenvalue of the L-shaped region. The first two curves report 
the difference between the values obtained with Richardson extrapolation of 124 — k grids, 
respectively using scheme iii and iv, and the precise value that we have obtained with the 
MPS; the last two curves report the difference between the values obtained with Richardson 
extrapolation of 124 — /c grids and the values obtained with Richardson extrapolation of 
124 — k — 1 grids, respectively using scheme iii and iv. This difference essentially provides 
the number of stable digits achieved. In the first four curves the FD eigenvalues are obtained 
with an accuracy of 220 digits; the last curve is analogous to the second one, limiting the 
accuracy of the FD eigenvalues to 60 digits 


60 

50 

40 

■ 1 ) 

■8 30 
o 

20 

10 

0 

Figure 3: Correct digits of the lowest eigenvalue of the L-shaped membrane obtained with 
Richardson extrapolation using a set of FD grids with a smallest spacing hmin = 

Notice that the number of grids used for a given hmin depends on hmin itself (see Fig. n. 
The dashed curve is the fit f{n) = 7.23166 -|-0.383229n — 20iilZ6 "pjjg pp eigenvalues used 
in the extrapolation were computing using 220 decimal digit floating point arithmetic. 
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Figure 4: Optimal number of FD grids used for a set of FD grids with smallest spacing 
hmin = The dashed curve is the fit g{n) = 0.542696n + 4.20652. The FD eigenval¬ 

ues used in the extrapolation were computed in multiple precision floating point arithmetic 
with a precision of 220 decimal digits. 



Figure 5: Error in the first eigenvalue of the L-shpaed domain using the diagonal Pade- 
Richardson Extrapolation p(^424) horizontal line corresponds to the minimal error 

obtained with the Richardson extrapolation, corresponding to the minimum in Eig. [21 
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this case, the fourth scheme is the appropriate one and the asymptotic expansion is obtained 


E{h) « 15.19725192645434327487838213300054590 06 + 3.18 • 10~^^ 

- 12.565568615260003775714180770 53 

- 2.2529040988480935561491817 46 - 9.9 • 10"^® 

+ 3.932508901213713526500 07 + 1.1289729308101123792 42 

+ 0.950164117523872693 21 - 1.36902891120799 30 /i® 

- 0.0740361191169 66 + 0.772459685647 50 + ... (40) 


Notice that in this case we have used the less precise set of FD values, which were 
computed only in 60 digit floating point arithmetic: the eigenvalue of the first excited state 
is now reproduced with “just” 37 correct digits. 

This result clearly shows that the coefficients of the terms and must vanish: in 
particular it is easy to understand the absence of since the mode that we are calculating 
is the fundamental eigenmode of the desymmetrized region obeying Dirichlet boundary 
conditions ony = x. In this case the reentrant corner is 7r/a = 37r/4 and therefore 2a = 8/3. 

With this simple observation, eliminating 4/3 and 10/3 from the exponents used in the 
extrapolation scheme, we are able to obtain 3 more digits of E 2 

|^(RE)_^(MPS)| = 40-8 

Even more digits can be obtained using the Fade-Richardson scheme, without the expo¬ 
nents 4/3 and 10/3: in this case 

I „(PRE) iTMPSjT ^ 45.8 

1^2 ^2 I 


4.2 H-shaped domain 

We now consider a domain with the shape of H, displayed in Fig. [31 originally studied by 
Donnelly m using the method of particular solutions (MPS) and finite differences (FD). 
As we have already mentioned in the previous section, the author conjectured that the FD 
eigenvalues, corresponding to a given grid spacing h, behave as 

Eih) = E{0) + ah^!'^ + bh^ + + ... (41) 

where E(0) is the corresponding eigenvalue of the Laplacian in the continuum and the 
exponent 4/3 is determined by the presence of a reentrant corner 37r/2 [TUI 121) . 

As for the L-shape, we want to obtain a precise estimate of the lowest eigenvalues for this 
problem, using a sequence of FD eigenvalues, obtained for different grids. Notice that the 
eigenfunctions of the Laplacian on this domain can be classified according to four different 
symmetry classes, even-even, even-odd, odd-even and odd-odd with respect to reflection 
about the x and y axes. By working separately on the modes belonging to each class, the 
computational complexity of the problem can be reduced and finer grids can be studied. 
Our present analysis, in particular, is limited to the even-even modes. The spacing of the 
grid is chosen so that the border of the H-shaped is sampled exactly and it corresponds to 

hk = 3/2/(9 -I- 3(fc — 1)), with fc = 1,2,_ We have calculated the first 25 eigenvalues of 

the even-even modes of the H-shape with a floating point precision of 60 digits, for the grids 
corresponding to k = 1, 2,..., 40. 
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Table 2: Lowest 25 eigenvalues of the L-shaped domain obtained with the MPS using 425 
points evenly spaced on the border. The eigenvalues marked with f are known exactly; the 
first eigenvalue, marked with =i<, has been obtained using the MPS with 545 points. 


n 

p(MPfc!) 

1* 

9.639723844021941052711459262364823156267289525821906456109579700564036 

2 

15.197251926454343274878382133000545900777179939609 

3+ 

27r2 

4 

29.521481114144883298220387998949268230835182037083 

5 

31.912635957137762200327505645485619891180683442197 

6 

41.474509890214922338810104064796906887679915692804 

7 

44.948487781351230152829670239630032397049780134665 

8+ 

57r^ 

9+ 

57r^ 

10 

56.709609887385120714216741638492259079610565870838 

11 

65.376535709845878509384400627738811907191161706097 

12 

71.057755648513529930798223378765313509589316160842 

13 

71.572679680336556014706999077329408038228565031443 

14 

87r2 

15 

89.301668351960185629207557215836143584908527108716 

16 

92.306906763049247832266397297040944898714305036279 

17 

97.380722646021860253461536778106579066564981169123 

18 

IOtt^ 

19 

IOtt^ 

20 

101.60529408377871548543481415097538087072356189211 

21 

112.36860922562569413546584663077376004912074741174 

22 

115.52017309466770886932756039014897616475657545671 

23 

137r^ 

24 

137r^ 

25 

130.11902885096790256577606801292831058988583848246 
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Table 3: Correct digits of the first 25 eigenvalues of the L-shaped domain, obtained applying 
the Richardson and Richardson-Pade extrapolations to FD eigenvalues. The values marked 
with the t correspond to eigenstates of the square. The first eigenvalue has been obtained 
extrapolating the FD eigenvalues of 124 grids, obtained with a floating point precision of 
220 digits. 


n 

scheme 

1 ^ ^ i 

1 ^ ^ i 

parity 

iOglO ^{RE)_^{MPS)| 

^•^610 ^{PRE)_^(MPS) 

1* 

iv 

54.5 

67.5 

even 

2 

iv 

40.8 

45.8 

odd 

3t 

i 

62.9 

73.1 

even 

4 

iv 

37.1 

45.8 

odd 

5 

iv 

35.9 

42.6 

even 

6 

iv 

35.1 

42.2 

even 

7 

iv 

36.7 

44.5 

odd 

St 

i 

60.6 

73.9 

odd 

91' 

i 

60.8 

73.8 

even 

10 

iv 

35.2 

41.9 

even 

11 

iv 

34.5 

42.6 

odd 

12 

iv 

34.8 

42.2 

even 

13 

iv 

34.4 

42.6 

odd 

14t 

i 

60.3 

73.2 

even 

15 

iv 

33.4 

41.3 

even 

16 

iv 

30.8 

39.9 

odd 

17 

iv 

30.6 

39.3 

odd 

18t 

i 

60.3 

74.0 

odd 

igt 

i 

59.5 

73.9 

even 

20 

iv 

33.0 

40.7 

even 

21 

iv 

32.6 

40.0 

even 

22 

iv 

33.7 

42.6 

odd 

231' 

i 

59.6 

73.6 

odd 

24t 

i 

59.7 

73.2 

even 

25 

iv 

33.3 

43.4 

odd 
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Our results for the lowest eigenvalue should be compared with those of Donnelly [10] 

^(Donnelly) ^ 7 jggQggg ( 42 ) 

and, more recently, of Betcke and Trefethen [3] 

= 7.7330888559 (43) 

In Fig. Owe report the error over the first eigenvalue of the H-shape. The first two curves 
report the difference between the values obtained with Richardson extrapolation of 40 — A: 
grids, respectively using scheme iii and iv, and the precise value of Betcke and Trefethen |3]. 
However, since the results of Ref. |3] are not sufficiently precise, it is convenient to estimate 
the error using the difference between the values obtained with Richardson extrapolation 
of 40 — fc grids and the values obtained with Richardson extrapolation of 40 — fc — 1 grids, 
respectively using scheme iii and iv. This difference essentially provides the number of stable 
digits achieved. Notice that the second curve rapidly reaches a plateau, for k < 34, signaling 
that in this range the extrapolated results are more precise than those of Ref. |3] . 

The figure clearly shows that the asymptotic behavior conjectured by Donnelly in Ref. [10] 
is not correct; our best estimate of the fundamental eigenvalue corresponds to the last curve 
in Fig. |3(i.e. scheme iv) for k = 18: 

El = 7.7330888559426190667 (44) 

where all the digits are believed to be correct. 

In table [4] we report the approximate values of the first 24 eigenvalues of the even-even 
modes of the H-shape obtained using Richardson extrapolation. It is particularly interesting 
to consider the value for the mode 24, which has the lowest precision. The coefficients of 
the asymptotic series obtained from the Richardson extrapolation are (underlined digits are 
expected to have converged) 

E{h) 194.7347257248 53 -k 1.2880 50/t'^/^ - 2861.993 46/t^ 

- 25.761/t^/^ -f 515.6h^°/^ + 14691.3/^"^ -k ... (45) 

The coefficients of this series, although determined with less precision than in the cases 
discussed earlier for the L-shape, clearly suggest the presence of a smaller radius of conver¬ 
gence, which drastically affects the accuracy of the calculation. 
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Table 4: Lowest 24 eigenvalues of even-even modes of the H-shaped domain obtained using 
Richardson extrapolation with set iv (the sets marked with f are eigenstates of the square 
and are extrapolated using set i). 


n 

7J, (Richardson) 

1 

7.7330888559426190667 

2 

14.30522996107150163018552 

3t 

19.73920880217871723766898199975230227062739 

4 

33.0048892952083545188 

5 

37.2054234400574157525 

6 

46.2961910861973723751 

7 

58.7501048292892847997 

8 

63.113298546574958190 

9 

67.43457224647486521 

10 

85.80372978847046992 

11 

92.12485042399187898 

12 

95.7615825533281487 

13+ 

98.696044010893586188344909998 

14 

112.42755013401679304 

15 

122.557976404091254965 

16 

133.5364354179283 

17 

139.4282184592822 

18 

142.4312241050896 

19 

150.543062476658690 

20 

164.339040164448839 

21 

171.85972578742946 

22+ 

177.652879219608455139020837997770 

23 

180.46602205029118 

24 

194.7347257248 
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Figure 7: Error over the first eigenvalue of the H-shaped region. The first two curves report 
the difference between the values obtained with Richardson extrapolation of 40 — fc grids, 
respectively using scheme iii and iv, and the precise value of Betcke and Trefethen [3] ; the last 
two curves report the difference between the values obtained with Richardson extrapolation 
of 40 — A: grids and the values obtained with Richardson extrapolation of 40 — fc — 1 grids, 
respectively using scheme iii and iv. This difference essentially provides the number of stable 
digits achieved. 


4.3 Isospectral domains 

Consider the domains of Fig. El It is known that these domains are isospectral, i.e. that the 
eigenvalues of the laplacian on one domain coincide with those on the second domain, as 
proved by Gordon, Webb and Wolpert [mill]. The numerical calculation of the eigenvalues 
of these regions has attracted large interest, using different techniques; for example, Wu, 
Sprung and Martorell [39] have used finite difference and mode matching to estimate the 
first 25 eigenvalues of these domains; the most precise results have been obtained by Driscoll 
in Ref. m and by Betcke and Trefethen [3]. The result that Betcke and Trefethen report 
for the eigenvalue of the fundamental mode 

El « 2.537943999798 

is slightly more precise than the value reported by Driscoll. Moreover, Sridhar and Ku- 
drolli [34] have performed an experiment with microwave cavities of the form of the domains 
of Fig. El verifying their isospectrality H. 

In this case, we have applied finite differences calculating the lowest eigenvalues of both 
domains for 30 grids; the grid spacing is chosen appropriately so that the border is sampled 
exactly 0. Remarkably, the matrices obtained with finite difference for the two domains are 
also isospectral. 

In Fig. [9]we report the error over the first eigenvalue of the isospectral domains, while in 
TableElwe report our best estimates for the lowest 25 eigenvalues, obtained using Richardson 

® Readers interested in the topic of isospectrality should refer to the recent review paper of Giraud and 
Thas US]. 

“^With respect to the case of the L-shape, here the domains do not have any symmetry and only specific 
grids sample the border; this explains the smaller number of grids which could be used. 
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Figure 8: Isospectral domains 


extrapolation, with the same exponents as for the L. For the lowest eigenvalue we gain 5 
digits with respect to the result of Betcke and Trefethen 

El = 2.53794399979862045 (46) 

Moreover, even our poorest result, for the 25th mode, has two extra digits with respect to 
the result of Driscoll. 

In light of these results, we stress that the finite difference method can provide very 
accurate results, despite the common prejudices. In the abstract of the paper of Driscoll, 
for example, we read: ’’Furthermore, standard numerical methods for computing the eigen¬ 
values, such as adaptive finite elements, are highly inefficient”. 

A second comment regards the work of Wu, Sprung and Martorell, who calculated the 
FD eigenvalues for these domains for 3 grids and then used Richardson extrapolation to 
obtain better estimates. Incorrectly, they assumed that the FD results vary quadratically 
with the grid spacing, a behavior which is appropriate only for the modes of the square 
(modes 9 and 21). 

4.4 Square domain with a 45°-crack 

The domain represented in Fig. [10] is particularly interesting, since it contains a reentrant 
angle 0 = iTTjA, which is larger than the angle of the L-shaped domain. Additionally, the 
domain has no symmetry and therefore the numerical calculation is more demanding than 
for the case of the L and H shapes. This problem has been originally studied by Blum and 
Rannacher [3| and more recently by Yuan and He [30]) where the bounds 

35.631515 <Ei< 35.631522 

have been obtained. The result Ei ss 35.617 was obtained in Ref. [4] applying Richardson 
extrapolation to finite elements. 

In Table |6| we report the numerical approximations to the lowest 5 eigenvalues of this 
domain, obtained using the MPS with 356 points. The digits reported in the table are 
expected to be correct; in particular for the lowest eigenvalue we have 

El « 35.63151951719172309520548614207765698409 (47) 

In Fig.[Tl]we show a contour plot of the hrst four modes of this domain, obtained using 
finite differences with a grid with spacing h = 1/120, corresponding to a total of 12331 grid 
points. The solid blue lines are the nodal lines, while the dashed green lines are level curves. 
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Figure 9: Error over the first eigenvalue of the isospectral regions. The hrst two curves 
report the difference between the values obtained with Richardson extrapolation of 30 — A: 
grids, respectively using scheme iii and iv, and the precise value of Betcke and Trefethen [3] 
{El K, 2.537943999798); the last two curves report the difference between the values obtained 
with Richardson extrapolation of 30 — A: grids and the values obtained with Richardson 
extrapolation of 30 — A: — 1 grids, respectively using scheme iii and iv. This difference 
essentially provides the number of stable digits achieved. 



Figure 10: Unit square with a 45°-crack 
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Table 5: Lowest 25 eigenvalues of the isospectral domains obtained using Richardson extrap¬ 
olation with set iv (the sets marked with f are eigenstates of the square and are extrapolated 
using set i). 


n 

7j,(llicharclson) 

1 

2.53794399979862045 

2 

3.65550971352441826 

3 

5.17555935622451540 

4 

6.53755744376443310 

5 

7.2480778625641275588 

6 

9.20929499840321242 

7 

10.59698569133316780 

8 

11.5413953955859566289 

9t 

12.33700550136169827354311374984518891914212 

10 

13.0536540557280658 

11 

14.313862464291008706 

12 

15.871302620009314 

13^ 

16.941751687972089 

14 

17.6651184368431201 

15 

18.9810673876525993 

16 

20.882395043282328 

17 

21.2480051773728 

18 

22.23285179297328 

19 

23.711297484824032 

20 

24.479234069273887 

2lt 

24.674011002723396547086227499690377838284 

22 

26.08024009965984 

23 

27.304018921125 

24 

28.175128581453 

25 

29.569772913239 
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Figure 11: Nodal lines of the first four excited inodes of the unit square with a 45°-crack 


While the fundamental mode is nodeless, the remaining three states have one or two nodal 
lines which start on the vertex of the reentrant corner, thus dividing the original domain in 
two or more domains. Looking at the figure we see that for the second state the resulting 
sub-domains have a reentrant angle 9 = Ttt/S, while for the third and fourth states the 
sub-domains have a reentrant angle 9 — 77r/12. The dashed straight lines in the plot are 
tangent to the nodal line in the vertex. 

As a result of this observation, we speculate that the asymptotic behavior of the finite 
difference eigenvalue may contain the exponents 8/7, 16/7 and 24/70. 

We have calculated the lowest eigenvalues for this domain using finite difference with 60 
grids; the Richardson and Richardson-Pade extrapolations of these results, with the appro¬ 
priate exponents in the asymptotic series, should allow one to obtain precise approximations 
to the eigenvalues of this domain, as for the case of the L. 

In this case we have extrapolated the finite difference results using a series of the form 

E[h) = E{ 0 ) + Cl + C 2 h'^ + C3 -I- C4 -I- C5 -|- ce h'^ 

+ cr + cs + C 9 /I®®/" + CIO + CIO + cn + C 12 

+ Ci3 /i® -I- Ci4 h®^/^ -I- Cio -I- C16 -I- Ci7 -I- C18 /l®®'^’^ 

+ CIO + C20 + C21 ^ /^136/7 ^ ^20/7 

+ C25 + C26 + ... (48) 

where the coefficients are chosen empirically and include the ones mentioned earlier. 

It is interesting to check the numerical values obtained for the coefficients of the series 
(l48)l . using the Richardson extrapolation of the FD results corresponding to the last 30 

®In the case of the L-shape, the reentrant corner is divided in two halves by the line y = x for the modes 
that are odd: in that case, the nodal line is exactly sampled by the grid and therefore the exponent 4/3 is 
absent, while the first rational exponent is 8/3. In the present case the nodal lines are not sampled by the 
grid. 
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Table 6: Lowest 5 eigenvalues of the unit square with a 45°-crack obtained with the MPS 
using 356 points evenly spaced on the border 


n 

p(MPS) 

1 

35.63151951719172309520548614207765698409 

2 

54.19310844424629197411978585647040768914 

3 

73.63330812560383459483828674566950026083 

4 

104.3280904734882128897772035674716112638 

5 

124.5914636064409738708659060017320376707 



Figure 12: Error over the first eigenvalue of the unit square with a 45°-crack. The asymptotic 
series of Eq. (H51) has been used. 


Table 7: Correct digits of the first 5 eigenvalues of the unit square with a 45°-crack, obtained 
by applying the Richardson and Richardson-Pade extrapolations to ED eigenvalues. 


n 

1 ~ w i 

1 ~ w i 

|£.^RE)_^(MPS) I 

-*^0810 j^^PRE)^^MPS)| 

1 

2”2.18 

25.37 

2 

23.65 

23.87 

3 

22.00 

23.92 

4 

21.04 

24.22 

5 

20.85 

23.01 
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grids, for the modes above: 

Ei{h) « 35.63151952 + 22.47641559 - 71.03523727 + 6.078713368 

- 78.46323288 - 8.840565052 + 63.35756993 + ... (49) 

E 2 {h) 54.19310844 - 2.87 x 10"^^ - 164.3992546 - 21.20457267 

+ 1.03 X 10”® - 1.44 X 10”^ + 212.7295338 h'^ + ... (50) 

Eaih) ?» 73.63330813 + 3.52 x 10"^^ h®/^ - 260.5413126 + 8.56 x 10"^^ 

- 3.15 X 10"® 91.25393089 /1^4/7 + 222.794824 /1^ + ... (51) 

Ei{h) PS 104.3280905 - 2.12 x lO"^® /i®/^ - 668.8593013 /i^ - 3.38 x IQ-^® /i^e/^ 

+ 9.07 X 10"^/i^^/^- 39.10703889 /1^4/7+ 1997.967306 /i4 + ... (52) 

E^{h) Ki 124.5914636 - 2.6 x lO"^® /i®/7 - 766.4031071 - 13.2187842 /i46/7 

+ 1.17 X 10"® /i^2/^ - 0.00001758167793 /i^4/7 1901.063425 + ... (53) 

Clearly one observes that depending on the mode chosen, some of the coefficients are 
consistent with a vanishing value: these observations are summarized in Table [51 where the 
leading rational coefficients and the corresponding reentrant angle are reported for each of 
the first 5 modes. 

Table 8: Leading rational exponents of the FD series for the first 5 modes of the square 
with a 45®-crack, and corresponding reentrant angles. 


n 

leading exponent 

dominant angle 

1 

8 

i2L 

9 

/e 

■h 

Q 

2\ 


A 

k 


c; 

/e 



_ 1 _ 

_a_ 


4.5 Square domain with two slits 


Consider the unit square with two 1/4 slits, represented in Fig. |T31 This example has been 
studied in Refs. I1I13]. In this case the re-entrant corner is 27r, thus the leading exponent in 
the FD series is ai = 1. Eliminating the pollution of this contribution, Blum and Rannacher 
were able to obtain Ei = 35.728 for their finest grid. 

Consistently with our previous assumptions, we conjecture that the FD series has the 
form 

CO 

F;W=co + ^c,/i4 (54) 

1=1 

which is the typical form used in Richardson extrapolation. In this case. Bender and Orszag 
provide in [5] a nice explicit formula for the coefficient cq (Eq.(8.1.I6) of pag. 375 of their 
book), which in our notation reads: 


Co 


N 


E 


E;("+'=)(n + /c)^(-l)'=+^ 
k'.{N -k)l 


(55) 
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Figure 13: Square domain with two slits 


Our numerical experiments with this domain consist of two sets: 

• a set which contains the numerical approximation to the lowest eigenvalue of the 
domain calculated to 220 digits of accuracy using the COM, for 36 grids with h = l/2n 
and n = 8,10,..., 80; 

• a set which contains the numerical approximation to the lowest 50 eigenvalues of the 
domain calculated to 60 digits arithmetic using the internal Mathematica command 
Eigenvalue for 20 grids with h = l/2n and n = 8,10,..., 46; 

In table [9] we report the approximate values of selected eigenvalues of this domain, ob¬ 
tained using Richardson and Pade-Richardson extrapolation. The eigenvalue of the funda¬ 
mental mode is obtained using the first set of FD results, whereas the remaining eigenvalues 
are obtained using the second set. The digits reported in the table are believed to be correct. 
The table omits the eigenmodes of the square, for which the convergence is much faster. 


Table 9: Selected eigenvalues of the square with two slits obtained using Richardson and 
Pade-Richardson extrapolation of the FD results 


n 

-t^n 

-t^n 

1 

28.131367480845754755206 

28.131367480845754755206268 

3 

70.65038470368 

70.65038470368488 

5 

99.846759253895 

99.8467592538950 

7 

130.483305932580 

130.4833059325804 

8 

153.39663535893 

153.3966353589373 

10 

196.598428600514 

196.5984286005142 

13 

218.04116455831 

218.0411645583168 

15 

268.2038796851519 

268.2038796851519 

16 

272.5993876495 

272.59938764953 

17 

280.750584654 

280.7505846542989 

20 

348.460286264284 

348.4602862642840 

50 

750.8475130 

750.847513086 
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Figure 14: Nodal lines of the 50*^ mode of the square domain with two slits. 


Of particular interest is the fiftieth mode, whose nodal lines are the solid lines displayed 
in Figs. [141 Looking at the left plot, we are tempted to assume that a nodal line partitions 
each of the 27r reentrant angles into three angles of 27r/3, which would imply that the 
corresponding FD series would now have rational exponents. A simple analysis of the FD 
results however shows that this mode is also described by the series in eq. (I54|) . This behavior 
is consistent with the information delivered by the right plot in Figs. HU that reveals that 
in fact the nodal line do not end in the reentrant corner. In other words, the study of the 
FD series for a given domain, can also provide information on the behavior of the nodal 
lines of the corresponding eigenmodes. 

5 Conclusions 

In this paper we have showed that it is possible to obtain precise estimates for the eigenvalues 
of the negative Laplacian over particular domains in the plane by performing a Richardson 
extrapolation or a rational (Pade)-Richardson extrapolation of the results obtained with 
finite differences, where the exponents of the series are related to the reentrant angles in 
the domain. The problem of determining the series describing the behavior of the finite 
difference results from first principles is difficult and it seems that a theoretical study is 
still lacking. The problem is both challenging and interesting for the applications of finite 
differences in Physics, Applied Mathematics and Engineering are as numerous as the stars in 
the Milky Way. Quoting Kuttler and Sigillito, pag. 178 of m,” the exact form of the first 
several terms in the asymptotic formula for specific regions where no boundary interpolation 
is required is a nice problem at about the level of a doctoral thesis.” The fact that, since 
1984 this problem has not been yet solved suggests an even higher level of difficulty. 

In this paper we have pursued the less ambitious goal of identifying the series (i.e. the 
exponents) empirically and we have obtained particularly encouraging results. In the case 
of the L-shaped domain, for instance, the extrapolation of the results obtained with finite 
differences leads to a determination of the first 68 digits of the lowest eigenvalue. 

The knowledge of the finite difference series for a given domain allows a precise determi¬ 
nation of the numerical values of the ei gen values of that domain, making the finite difference 
method a powerful computational toollj. 

Here we stress the most relevant observations obtained from a careful analysis of the 
numerical results for the examples considered in this paper: 

®In all the examples that we have treated in this paper, we have been able to improve published results. 
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• The FD series appears to be an asymptotic series, as suggested by the particular 
behavior of the error; this does not limit the accuracy of the extrapolated results, if 
the largest spacing of the set is appropriately decreased, as more and more terms are 
added; 

• The example of the square with a 45° crack tells us that when a nodal line terminates 
in a reentrant corner, the corresponding FD series have exponents corresponding to 
the fractions of reentrant angles, even if the nodal line is not completely sampled by 
the grid (it is the behavior infinitesimally close to the corner that matters); 

• It is reasonable to assume that, for a given domain, the FD series corresponding to the 
different modes all are described by the same series (although for some modes some 
exponents could be missing for symmetry reasons - this is the case of the modes of the 
L which are also eigenmodes of the square, for which all the coefficients of all rational 
exponents vanish ); 

• If the observation above is correct, this means that one cannot have nodal lines parti¬ 
tioning the reentrant corner if the new exponent generated is not of the type already 
contained in the series! The case of the fiftieth mode of the square with two slits 
illustrates this behavior: the nodal lines stretch almost completely to the reentrant 
corner, although they do not join it! 

• We conjecture that the nature of the reentrant corners fully determines the exponents 
of the FD series and therefore different domains, containing the same reentrant angles 
should all have the same exponents (see for example the case of the L, of the H 
and of the isospectral domains considered in this paper); this makes Richardson (and 
Richardson-Pade) extrapolation practical even for complicated domains where the use 
of MPS can be problematic; 

• For the case of the L-shape and of the square with a 45° crack, our results also provide 
an independent check/validation of the corresponding results obtained using MPS; 
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